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^ Abstract 

We study a discrete time queueing system where deterministic arrivals 
have i.i.d. exponential delays The standard deviation a of the delay 
is finite, but its value is much larger than the deterministic unit service 
time. We find the bivariate generating function for the system, and we 
solve the resulting boundary value problem in terms of a power series 
1— 1 expansion in a parameter related to a^ 1 . We also prove the analyticity 

0^ of the generating function with respect to this parameter. The model, 

motivated by air and railway traffic, has been proposed by Kendall 1 36 [ 
(—1 and others many decades ago, but no solution has been found so far. 

B 



INTRODUCTION 



In this paper we study a one server queueing system defined in the following 
way. The z'-th customer arrives to the system at time 



> 
ON 

ON *i = 1 + & (i-i) 

on 

y—i where are i.i.d. exponential random variables, with density 

^ /fU |0 otherwise 

k*" In the limit A — » this arrival process weakly converges to a Poisson process 

k> whereas for A small but fixed the arrivals are negatively autocorrelated, 

$_j see [30] and references therein. In this paper we study the system for fixed A. 

In order to ensure tightness for the queueing process we assume throughout 
the paper an independent thinning to the arrival process: each customer can 
be removed/ deleted before joining the queue with independent probability 
1 — p. In this case it is not difficult to show that the arrival process converges 
weakly to the Poisson process as above, see again O30I. After Kendall (see 
below for his exact words) we will name this arrival process exponentially 
delayed arrivals. 

Service can be delivered by the unique server only at discrete times. The 
length of the queue nt at time t is the number of customers waiting to be 
served, including the customer that will be served precisely at time t, if any. 
Due to the thinning procedure, it is immediate to see that the traffic intensity 
of the system is given by p; see [30] for details. 

This model is motivated by the description of public and private transporta- 
tion systems including buses, trains, aircraftsfgj[30j|3i}[33] and vessels II28II34I , 
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appointment scheduling in outpatient services O [i6l [42J [43I crane handling 
in dock operations 820! [23] , and in general any system where scheduled ar- 
rivals are intrinsically subject to random variations. Preliminary results show 
that the model described above fits very well with actual data of inbound air 
traffic over a large hub, see [15] . 

This model is an example of a queueing system with correlated arrivals, a 
subject broadly studied in past years. There are many ways to impose a 
correlation to the arrival process. For instance, the parameters of the process 
may depend on their past realization, as in II2TI , or on some on/ off sources, as 
in B57J. Among queueing systems with correlated arrivals Markov modulated 
queueing systems (MMQS) are quite relevant. In MMQSs the parameters are 
driven by an independent external Markovian process, see jH [7} [19] [41} [49} [51] 
and references therein. As we shall see in Section [2] our model shares with 
MMQSs the property that one can define an external and independent 
Markovian process that drives the arrival rates. However, the output of 
this external drive also determines the evolution of the queue length. More 
precisely, our system can be seen as a single-server queue with deterministic 
service time and arrivals given by the reneging from an auxiliary queue, 
namely the customers that are late at time t - see below. Due to 
the memoryless property of the exponential delays, each customer late at 
time t may arrive in the unit time slot (t, t+i] indipendently and with 
the same probability q = e~ A . This means that the aforesaid renenging 
only happens at integer times and that customers perform synchronous 
independent abandonments, leading to binomial transitions in the number 
of late customers. 

In Section [2] we will see that our model can be described as a bidimensional 
Markov chain, whose components are the queue length and the number 
of late customers, respectively. There exist an extensive literature about 
two-dimensional Markov models. Many methods for attacking the problem 
are available under the assumtions of either a spatial homogeneity property 
and that one of the two marginal chains is finite, see [Hoi [27} [29) [40II44] [46] pfS) . 
Unfortunately, the Markov chain defined in Section [2] does not satisfy any of 
the mentioned hypothesis. 

When both components of the Markov chain are infinite but space homo- 
geneity is ensured, the problem is typically attacked by reduction to a 
Riemann-Hilbert boundary value problem. Boundary value problems repre- 
sent a broadly studied subject and several techniques have been developed 
in the last decades to solve them: uniformization technique [37!, conformal 
mapping Ji7}[i8ll26ll , compensation method [3] and power series approxi- 
mation (PSA) riraiT2"l|3"2")|3"8) . Usually, PSA is used to obtain the generating 
function in terms of a power serie in the load p, though different parameters 
may be used, as in [55] . A power series expansion in a parameter different 
from p is also the strategy we choose in this paper. 

The second constraint that our model fails to meet, namely the spatial 
homogeneity, is due to the reneging with binomial transitions mentioned 
above. This kind of transitions are often encountered in Mathematical Biol- 
ogy [6[ [14J [22] . The functional equation for the boundary value problem | |2.8) 
below and its counterpart in P3] [24) [33] are radically different due to the 
mismatched dimensionalities of the systems, yet it is still possible to mark 
some analogies. The most important is that in both equations the right hand 
side exhibits the generating function computed in a convex combination in 
the parameter q, that is the probability of each independent abandonment; 
this feature will play an important role in the proof of Theorem [2] Other 
examples of chains with binomial transitions may be found in El |5) |47l [58| - 
Although some features of the model may be found in recent literature the 
problem studied in this paper was proposed and intensively studied for 
its remarkable importance by the founders of queueing theory in the '60s. 
The appearance of the stochastic point process can be traced back to 
Winsten's seminal paper fjp?) . Winsten named such a queueing model late 
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process and obtained results for the special case where customers can never 
be delayed more than two time-units, i.e. G [0,2]. At the end of J36) , a 
discussion by Lindley, Wishart, Foster and Takacs on Winsten's results is 
attached. According to their words, Winsten's paper has to be considered the 
first treatment of a queueing model with correlated arrivals [56] pg. 22-28]. It 
is worthy to note that in the framework depicted above Winsten surprisingly 
got a geometric distribution for the stationary queueing length. During the 
discussion of Winsten's paper, Wishart [56] pg. 24] states the following: 'If we 
suppose that the £'s have the exponential distribution, then I think that we 
will no longer have the geometric distribution: but what we will have I do 
not know'. In this paper we give a contribution representing a step towards 
to the answer of that question. 

The same problem was investigated also by Kendall in [36 pg. 12], where 
he remarked the great importance of systems with arrivals like | |i.i) : '[...] 
perhaps too much attention has been paid to rather uninteresting variations 
on the fundamental Poisson stream. As soon as one considers variations 
dictated by the exigencies of the real world, rather than by the pursuit of 
mathematical elegance, severe difficulties are encountered; this is particularly 
well illustrated by the notoriously difficult problem of late arrivals.' Kendall 
also provided the following elegant interpretation: if the random variables 
are non negative, then the process defined in is the output of the sta- 
tionary D/G/oo queueing system; in particular, if the £,'s are exponentially 
distributed then the process defined in is the output of the stationary 
D/M/oo queueing system. 

Some years later, Nelsen and Williams obtained in J43H the exact characteri- 
zation for the distribution of the inter-arrival time intervals and correlation 
coefficient between successive inter-arrival time intervals when random vari- 
ables £j are non-negative. They also gave an explicit expression of these 
quantities in particular case where the £,'s are exponentially distributed. 
After that only approximations of e.g. JT3] [32] , and numerical studies 
of its output, see for example B4I Ig] [50) , seem to have appeared in literature. 
In particular, in [30] we studied in a self-contained way an arrival process 
defined as in {l.l) , assuming for a compact support distribution. We also 
proposed an approximation scheme that keeps the correlation of the arrivals 
and is able to compute in a quite accurate way the quantitative features 
of a queueing system with this kind of arrivals. Thus, to the best of our 
knowledge, a queueing system with arrivals described by still remains 
an open problem and the best results obtained so far are due to Winsten. 
In this work we present an equation for the bivariate generating function of 
the probability distribution of the queue length of the system in the case of 
exponentially delayed arrivals, and we write a formula leading to the explicit 
expression of the generating function as of a power series in the parameter q. 
The paper is organized as follows. In Section [2] we derive step by step an 
equation for the bivariate generating function. In Section [3] we develop 
an iterative method to compute the explicit expression for the generating 
function in terms of a power series. In addition, we discuss some aspects of 
the solution and presents in detail the first terms of its expansion. Section I4] 
contains a proof of the analiticity of our generating function with respect 
to the parameter q. Section [5] is devoted to the discussion of some open 
questions and some conjectures on the system. We also compare the queue 
length distribution obtained using the truncated power expansion with an 
empirical distribution obtained via simulations of (l.l) . 

2 CONSTRUCTION OF THE BIVARIATE GENERATING FUNCTION 

The process fif, describing the length of the queue at time f, is governed by 
the relation 

n t+1 = n t + m ( f (f+1 ] - (1 - S ntfi ) (2.1) 
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where m(t,t+i] is the number of arrivals in the interval (t, t + 1], and the term 
1 ~~ dn t ,0 represents the fact that if at time t there are clients in the queue then 
the first of the queue is served, decreasing the length of the queue by 1. 
It is evident that the quantity m^ ft+1 j depends in general on the whole 
previous history of the system. If for some large value of T we have that 
?W(- S;S+1 ] = for any s E {t — T, t — T + 1, t — 1}, then we will have a great 
probability that m^ t f+1 j is large; conversely if in the recent past the values 
of WZ( S s+1 j have been large, then m^t+x] is expected to be smaller. This is 
equivalent to claim that the arrival process is negatively autocorrelated, for a 
proof of this claim see [30]. 

It is then understood that the evolution ( |2.i) does not depend only on the 
present value of nt, indeed the memory of the process is infinite since T can 
be arbitrarily large. Nonetheless the delays are exponential, i.e. memoryless. 
Let us denote with If the number of customers that have not arrived yet at 
time t, that is to say 

k=\{0<i<t s.t. £>f-z}| (2.2) 

Next, define p = f^(t)dt = 1 — e~ A and q = e~ A . Then, given the value 
of U, the random variable m^ t >t+1 i is binomially distributed. More precisely, 
using the notation m^ f f+1 i = nit for the sake of simplicity, we have that 

P(m t =j\U = l)= ([)p>1 H = hl < 2 -3) 
if the customer at time t has been deleted by the thinning procedure, while 

p(m t = j\u = i) = C^ 1 } p'i 1+1 -' = Hw (24) 

otherwise. So we have 

P(m t =j\U = l)= Q pV-j (1 - p) + (' + ! ) P V +M P (2-5) 

= b j,l (1 — p) + bj,i+i p (2.6) 
Hence, let us assume that the state of the system is determined by If and by 



the length of the queue nt. Relations 1 2.3 1 and 1 2.4 1 ensure the fact that the 



number of arrivals does not depend on the previous history of the system, 
because nit depends only on It . The evolution of the pair (nt, h) is therefore 
Markovian. The resulting Markov chain is clearly ergodic, therefore the 
stationary measure P n \ exists and it is unique, provided p < 1. 
Consider the following bivariate generating function 

P(z,y) = £ P n , /Z y \z\,\y\ <1 (2.7) 

n,l>0 

We are now ready to prove the main result of this section: 

theorem 1 The bivariate generating function P(z,y) defined in satisfies 



P(z,y) = P[yq + ZV ] + 1 9 [(z - l)P(0,yq + zp) + P(z,ya + zp)} (2.8) 



Proof. We can write a set of equation for the stationary probabilities P n \ 
computing the probabilities to have Pn t+1 ,l t+1 m terms of the various P nt ,l t , 
and then neglecting the time dependency. We obtain 

Po,i = p(Ptj-i + Po,i-i)h,i + (1 - p)(Pi,i + Po,i)\i (2.9) 
Pi,i = P [{P\,i + Pq,i)Ki+\ + Pz,i-ih,i] 

+ (1 - P) \{Pu+\ + Po,i+i)Hi+\ + Pi,i\i\ (2.10) 

Pi,i = P [(Py+i + Po,l+l)h,l+2 + P2,M,l + P3,l-lh,l] 

+ (l - 1°) [{P\,l+2 + Po,l+2)h,l+2 + P2,i+th,m + P3,l\l] (2-11) 
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where bjj are defined in fe.j) and ^2.4) . Defining 

Pl(z) = E F n,lZ n 
n=0 

we can rewrite the system above in the compact form 
z — 1 

P( z ) = h,l + ZP 0< 1 + Z 2 P ,/+1 &2,/+2 + ■■■■] 

+ ^[P/_i(z) b 0J +zPi(z) b v+1 + z 2 P l+1 {z) b 2/l+2 + ....} 
z — 1 

+ C 1 - P) — =— [ p o,i Kl + zPo,i+i h,M + z2p o,i+2 h,i+2 + ■■■■] 

+ 1 ^P[P 1 (Z) b 0/l + zP /+1 (z) &y +1 + Z 2 P ;+2 (z) b 2/l+2 + ....} (2.12) 

Eventually, defining bj(z) = X^T=o ^n,l z " we arrive to b\{z) = +zp) 1 . Using 
the latter expression into | |2.i2[ | we get □ 

This equation represents a non trivial boundary value problem that will 
be studied in the next section. We conclude the present section with the 
following 

Remark 1. The marginal distribution for /, i.e., P(l,y) can be obtained by a 
recursive substitution of P(l,y) in the expression < |2.8} . We obtain 

00 

P{i,y) = Y\[i + P q k+1 {y-i)\ (2.13) 

/c=0 

This infinite product is well known in combinatorics with the name of 
q-Pochhammer symbol of the pair (p(y — 1), q). It is the ^-analog of the 
descending factorial (a.k.a. Pochhammer symbol), and when p(y — 1) is 
substituted by 1 it represents the well known Euler's function. It is worthy to 
outline that in the circle | q \ < 1 the ^-Pochhammer symbol is analytic for all 
the values of y. 

3 EXPANSION IN POWERS OF q 

We look for the solution of the boundary value problem in the class of 
the functions which are analytic in the region 

{{z,y,q) G C 2 x [0,<p] : \z\ < 1, \y\ < 1, < q < <p) 

If a solution is found in this class then by the uniqueness of the stationary dis- 



tribution P„ ; it must be the unique solution of 1 2.8 1. Consider the following 
power expansion 

P(z,y) = £^pW(z,y) (3-i) 

it>0 

which is convergent for < q < (p. An estimate of (p is given in Section^ The 
function PW(z,y) = [^]P(z,y) is the A:-th coefficient of the power expansion, 
namely 

Let v = z + q{y — z) then the equation \2.S) for the bivariate generating 
function can be rewritten as follows 



P(z,y) = (1-p + pv) 



mu) + PM-P(M 

z 



(3-2) 



This can be combined with a Taylor's expansion of P(z, y) around y = z 

P(z,y) = E^^|y>(z,y) (3-3) 
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For brevity of notation we will denote with the symbol 3yP(-,yo) the y'-th 
order derivative taken with respect to y and evaluated at y = yo- From 
we get the following expression for the /-th derivative of P{z,y) 



dyl 



P(z,y) 



y=z 



d' y P{0,: 



±P{z,z)-±P(Q,z) 



(l + p(z-l)) 



+ jpq J 



d ! y 1 P (z,z)- d[ J l P(0,z) 
z 



(34) 



Some remarks are now due. 

Remark 2. If q = then the number It of customers not arrived yet at time t 
is identically zero; thus we can infer that the coefficient of the zero-th order, 
PW(z,y), is independent of y. 
Remark 3. For y = z equation becomes 

1 - p + pz 



P(z,z) = P(0,z) 

from which it follows 

P(0,z) 

Let us define the functions 



1-p 

P(z,z) -P(0,z) _ P(0,z) 



4(z) = |- (P«(0,,/H 



: 1-p 
pW(z,y)-P«(0,y)' 



(3-5) 



(3-6) 



(3-7) 



for any j, k > 0, with the agreement that j = means no differentiation and 
that fl^ (z) = whenever j < 0. Substituting into and rearranging 
in powers of o we get 



An Z 



pW(0,; 



p(' c )(z,z) -P«(0,z) _ P( fc )(0,z) 



z 1 — p 

Remark 4. For any fl > the boundary conditions are 

Jp(0, 1) = 1 - p LzffZe's Law 
1 P(l, 1) = 1 Normalization 

From fe.i) we have that 

l-p = p(0,l) = ^/pW(0,l) V0<o<l 

/c>0 

When q = equation ( |3.io[ ) becomes 

P( 0) (0,1) = P(0,1) = l-p 
and from Remark [2] we get 

P(°)(0,y) = l-p Vy 
Equations fe.io) and I j3.11) then yield to 

PW(0,1)=0 VA:>1 



V Jfc > (3.8) 



(3-9) 
(3.10) 

(3-«) 
(3.12) 

(3-i3) 



We are now ready to write the fundamental theorem of this section that will 
allow us to determine the coefficients P( k \z,y) of the power expansion ( fe.i) . 
theorem 2 The coefficients P^ k \z,y) satisfy 



P(°)(z,y) = l + p(z-l) 



pW(z,y) 



7=1 

+ (l + p(z-l)) P *j°' z) Vfc>l 



(3-i4) 



(3-i5) 
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Proof. Equating 1 3.1 1 and 1 3.3 1, and using 1 13.4) we get 

d'p{z,z)-d' v p{o, z y 



E^ w (z,y) = E 

k>0 j>0 



[1+P(Z-1)) \^yP(0 /Z ) + 

1 



( i-l , ^ P(Z,Z) -d ] y P(0,Z 



7! 



(3.16) 



Using the analyticity we can now substitute fe.i) into the right-hand side 
of 1 3.16) and rearrange the expression to find that 



ir>n ;>n /• t>n 

(3-17) 



k>0 j>0 >' i>0 

Equating the zero-th order terms and using (I3T8J we find 



P<»), z , y )Hi + K--D]^ 



(318) 



and ( 3.14) follows from ( 3.12) 

Consider now k > 1 and operate the change of indices i + j — k in 1 3.17) 



k>0 k>0 ;=0 ' - 

(3-19) 

/c 



= E*< 



t>0 



E {/?£>(z) + [l +p(z - l)]^(z)} 



+ ( 1 + P(2 -i))^ 



(3.20) 



where from [ 3.19) to [ 3.20) we have used fe.8\ and the fact that fly(z) is 
identically whenever j < 0. Equation 1 3.20) easily gives □ 



Remark 5. Line 1 13.15) gives a recursive relation for PW(z,y) involving the 
derivatives of P^(z,y) only for i < fc and PW(0,z). In the following we 
show how to find an explicit expression for p(' c ' (z, y) using I 3.13) and I 3.15) . 

Remark 6. Note that P'") (1, 1) = 1. The second line of ([3.9) then infers 



PW(1,1)=0 Vfc>l 
Remark 7. From (3.7) , ( |3 .8| > and we easily see that 

«°(z) = V 
where m is the usual Kronecker 's delta. 



(3.21) 



(3.22) 



3.1 General form of P"(0,z) 

Let us define the following family of functions 



4( Z ) 



;p£>(z) + [l + p(z-l)]«*->(z) 

[l+p(z-l)] 




According to 1 3.23) we can rewrite line ( 3.15) as 



pW(z,y) = E^ Z i 2); Ak 
j=i h 



if 1 < j < k 
if ; = jfc = 
otherwise 



a;(z) + (i +/) (z-i))^M 



(3-23) 



(3-24) 
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Then we have the following result 
lemma 3 Foranyk>l, 

P«(0,2 



j=l J' 



(3-25) 



Proof. From 1 3. 13 1 and 1 3.24) 







which yields to 



pW (0/ o) = -£jii 



(3.26) 



7=1 



The value of pW(0,0) can now be used to computed PW (0,z). From ([3T24J 
again, we obtain 

P«(0,z) = £^4(0) + P«(0,0) 

;=1 7- ;=1 7- 

□ 



Combining equations ( [3.24) and I 3.25) we get the general form of pW (z,y) 



P«(z,y) = £ 
7=1 



(y- z > „* 



_ , , 1 + p(z - 1) — 1 . k/n . 



Remark 8. By means of fe.22) we have that, for k > 2 

A\{z) =p4_i{z) + [l+p(z - 1)] a&z) = 



(3.28) 



(3- 2 9) 



Therefore for any k > 2 we have that in formulas 1 3.15[ |, [ [3.24) and 1 3.25) the 
index j actually runs from 1 to k — 1 . 



3.2 First order in q 



For k — 1 equation 1 3.28 1 gives 



u , , l+p(z-l) 



p( 1 )(z,y) = (y-z)Al(z) + ^ 



(z-l)A{(0) (3.30) 



where Aj (z) = p by virtue of 1 3.22) . 
Thus we find 



P( 1 )(z,y) = (y-z)p+^ 



l+p(z-l) 



p(z-l) 



(3-3 1 ) 



Note that PW (1, 1) = according to I3.21 1. Differentiating with respect to y 
we find 



a\{z) =p 

«J(z)=0 V;>2 



Using 1 3.23) we have 



(Z)=p(l+p(2-l) i 



2-p 



(3-32) 

(3-33) 
(3-34) 



(3-35) 
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3.3 Second order in q 

From 1 3.24) , 1 13.25) and Remark [8] we get 



p( 2 )(z,y) = {y - z )A\{z)+ l+ P { \ l \ z-l)A\{Q) 



(3-36) 



which gives 



p( 2 )(z,y) = (y-z)p 



l + p(z-l) T 
1 ±? (Z f 1) p( Z -l) 



(3-37) 



Note again that P^ (1, 1) = according to 1 3. 21) . Differentiating we obtain 



ag(z) 



1 



1-p 1 



1) 



,2/ 



(3-38) 

(3-39) 
(340) 



1 + 10(2-1) 2 



aj(z)=0 V;>2 
Using 1 13.23) we have 

A?(z)=/ 
A3(z)=2p 2 

Remark 9. Since pW (z, y) and p( 2 ) (z, y) are linear in y, we have that 

At (z) = 4(z) = • • • = Aj +1 (z) =0 V; > 3 (3.43) 

= 



(3-4i) 
(342) 



^(z)^(z)^...^A; +2 (z) 



V; > 3 



(3-44) 



In particular, we see that the lowest k such that P( k '(z,y) is more than 
quadratic in y is k = 6. Such a value could be directly deduced from the 
memoryless property of the delays. The most likely situation leading to 
It = 3 (see | |2.2) above) is that none of the customers expected to arrive at 
times t — 1, t — 2 and t — 3 has arrived yet. This event has probability of 

order q 6 . Similarly the lowest k such that p( k '(z,y) is of order m in y is 

^ _ m (m+l) 



3.4 Third order in q 

From fe.24^ [ 3.25) and Remark [8] we get 

A\{z) + {y-z)A{{z) 



P {3 \z,y) = ^ 2)2 > ■"''•< 



l+p(z-l) 



l-p 

Thus 

pW(z,y) = (y-z) 2 p 2 + (y-z)p 
, l+|0(z-l 



[(z-l)A?(0) 
l+p(z-l) (2 



A1(0)] (345) 



(1-P) 2 



1 



V(z-l) 



•P 2 



g 

(1-P) 2 



+ p(z + l) 



(346) 



Note again that P^ (1,1) =0 according to < 3-2i| ). Differentiating we obtain 



«§(z) = I ^( z2 - 1 ) + r 
«f(z) = 2p 2 (z-l)+p 

fl3( Z )=2p 2 

«J(z)=0 V/>3 



1) 



l-p 2 



P 

(1-P) 2 



(347) 

(348) 
(349) 
(3-50) 
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Using 1 3,23) we have 



A\{z)=p 
A 4 2 (z)=2p 



l+p(z-l) 3 + 2p(z-l) 



l+p(z-l) 



(3-5i) 
(3-52) 



4 ESTIMATE OF THE RADIUS OF ANAL YTI CITY (f) 



In this section we provide an estimate of the radius of analyticity of fe.i) . 
Combining the definition ( [3.7) of the functions fly (z) with formulas 1 13.23) 
and 1 3.28) we easily derive the following recursive relation, valid for I > 1 



7)- 



-Z-l, 



z-1) 



/=/+! V l > 



;>«^(0) + (l-p)^(O) 



Mi 



( , , flf-'(z) — a*-' (0) \ , , 
+ IP «]--[ (0) + ~ + paf (z) 



z 

fc __«*-<(z) -«*-<(())' 



(4-i) 



while for / = through l fe.8) , [ 3,23) and I3.25 1 we get 



«g(z) = E 



1 z' -1 



ipa k r>(o) + (i-pyr>(o) 



pip- 

lemma 4 For flZZ A:, Z,wz e N flnd |z| < 1 



(4-2) 



dz m ' v ; 



(4-3) 



a;/zere C = max je 2 , jz^ }■ 



Proo/ By induction. For k = 0, 1 1 14.3) is clearly satisfied. In addition, from 
Remark [8] it suffices to prove I I4.3) for / < A:. Fhus let us assume that ( I4.3) is 
fulfilled in the polyhedral lattice 

P fc = {(/,/,«) e N 3 , /' < /c, Z < /'} 
and show that I I4.3) holds in 

P fc+1 =F jt U {(/,/, m) e N 3 ,; = fc + 1, / <fc + l} 



First note that from ( 3.23) and I3.28 1 it follows very easily that a\(z) is 
a polynomial in z with degree not larger than k. Fhus, using the inductive 
hypothesis and a Faylor expansion we have that for V m < k 



d m flf(z) - flf(0) fP 



dz" 



dz m 

~ k 



dz! 



7<(0) 



z j-m-\Mk(Q) 

dz> ' w 



E 

L;'=i 

<C k (k-m-l 



(44) 



and again m = means that no differentiation is taken. 



10 



From §4.2) we now have 



2) 



;! (1 - 


P) 


2ic k 


-i 


j\ (1 - 


P) 


lei 





(jpC k -' + (l-p)C> 



< c k — — 

- C(l-p) 



(4-5) 



If C > ^- then (£5) infers |flg(z)| < C*. 
From ^4.1) instead we have 

< t jj^jyiC k -i + 2lpC k - l + pC k - 1 +2(1 -p)C 



a\(z) 



•k-l 



i=l+l 



< (21+2- p)C k ~ l + 
2-p + 2eh (l+ I 



lc k - 1 V l -M 



< c 



c* 



(4.6) 



and I I4.5) infers \a k (z)| < C k provided 1 
proof of I I4.3) for m > 1 is similar. 

< = c-\ 



The proof .... ... 

Lemma [4] implies then that 



C > e 2 . 
the analyticity radius is <p 



□ 



5 OPEN PROBLEMS AND NUMERICAL APPROXIMATIONS 



Although equation I 3.15) represents the complete solution of the problem, 
its use is quite involved, due to the fact that it is a recursive relation on 
functions. It would be nice to have an explicit expression for PW(z,y), 
or equivalently PW(0,y), involving only a (linear) recursion on a set of 
constants. Alternatively, it would be nice to have a similar tool to compute 
relevant quantities like the mean value of the queue and its higher momenta. 
Both questions remain unanswered, but it is clear that the problem has a 
deep combinatory structure. Here we want briefly mention that the evidence 
of such combinatorial structure emerges also in the easy case of the marginal 
distribution represented by P(l,y). According to (3.15 1 we can write for 
pW(l,y) the following expression 

p {k) (i,y) = E il ^P L frvld) +<f (5-!) 

;>0 J- 



And we have that 



«;.(!) = ^P«(i,y) 



P-P'D(lj) 



(5-2) 



where D(l,j) is the number of partition of the integer / in terms of j distinct 
parts, i.e. 

D(Lj)= E 1 (5-3) 
h<...<ij 



Hence ( fc.i) can be rewritten as 

P«(l,y) = E (P(!/ " 1))y 



D(k,j) 



j>0 h 

which proves that is equivalent to 

D(k,j) = D(k-j,j-l) + D(k-j,j) 



(54) 



(5-5) 
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Last equality can be easily verified by means of Young diagrams, see for 
example I54I. 

Finally, we want to state a conjecture that seems to be very reasonable looking 
at the first terms of the expansion in q. We conjecture that the general form 
of P{z,y) should be 

P{z,y) = J! ( 1 + pq k+ \y -z)+q k Y: h,^ " *) ) (5-6) 

*:>0 \ ]>\ ) 

with |8fc 1 constants to be determined. 

Another aspect of this model worthy of being studied is the convergence 
to the stationary measure of the system. It is very likely that the system 
presents the phenomenon of the cutoff ( for a discussion of the cutoff see 
for instance [[39] [33] and references therein). A proof of such cutoff can be 
probably achieved with the ideas in |39j|, using the parameter oct = nt + It- 
When the system starts from a very high value of tit, it is easy to prove (see 
[30] for more details) that at decreases by one when the customer t is deleted 
by the thinning procedure, and it remains constant otherwise. Moreover in 
this case «f is very far from its stationary distribution, leading to a large total 
variation distance of the distribution with respect to the stationary one. Once 
oct arrives to the values on which its stationary distribution is supported, the 
convergence to the equilibrium should be very fast. This will be the subject 
of further studies. 



The last part of this section is devoted to some numerical aspects of the 
model. In Section kl we explicitly found the first four coefficients p( k >(z,y) of 
the expansion fe.i} . A question that may naturally arise now is 'how good is 
approximation found this way?'. It turns out that the answer depends on the 
values of p and q. 

From 1 3.31) , 1 13.37) and 1 3.46) we see that the convergence of the truncated 
series YJk=o ( z '¥) 9* ^° -P( z '!/) ma y ^ e slow when p is close to 1 due to the 
presence of factors (1 — p) -1 and (1 — p)~ 2 - On the other hand we obviously 
have that the smaller is q the faster the aforesaid convergence. 
Numerical simulations can give an idea of the quality of the approximation 
of the truncated expansion. The probability distribution of the queue length 
can be found as the marginal of the empirical distribution P n \ or by means 
of the following formula 



1 8" 



(z,y) = (0,l) (5V) 



n\ dz n 
1 N ?)" 

I E^4 p( *^y),.„ wn +o( q N ) (5.8) 



„l L-i I far, 
k=0 °~ 



(z,y)=(0,l) 



Figure [5] shows a comparison between the distribution of the queue length 
empirically obtained via numerical simulations of an Exponentially Delayed 
Arrivals process and the distribution obtained by means of < |5.8) . 
It is immediate to see that if the traffic index p is close to 1 then even 
small values of q are sufficient for huge discrepancies to arise between the 
simulated and the theoretical distributions. Conversely, if the average load of 
the system is kept at a moderate level then we see that the model is capable to 
deliver a very good accordance with the simulated data for scenarios where 
each customer has a great probability to arrive out of its pre-scheduled slot. 
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Figure 1: Comparison of P n = XjP„/ for interesting scenarios. For each 
different couple of values (p, q) the 'Empirical' serie represents the empirical 
distribution while 'Theoretical' is obtained by means of l |5.8) . 
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